#include <iostream>
#include <cmath>
#include <string>
#include <fstream>
#include <sstream>

#include <TROOT.h>
#include <TFile.h>
#include <TTree.h>
#include <TRandom.h>
#include <TMatrixD.h>
#include <TLorentzVector.h>
#include <TH1D.h>
#include <TH1I.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TGraphErrors.h>

#include "GeomSvc.h"
#include "SRawEvent.h"

using namespace std;

int main(int argc, char *argv[])
{
  GeomSvc *p_geomSvc = GeomSvc::instance();
  p_geomSvc->init(GEOMETRY_VERSION);

  SRawEvent *event = new SRawEvent();
  int nSeeds, nSeedsX, nSeedsY;
  double ax[50], ay[50], bx[50], by[50];
  int xIndex[50], yIndex[50];

  double time;
  int nCandidates;
  double chisq;
  int charge;
  int nRecHits;
  int iRecHits[50];
  int nHodoHits[3];
  int detectorID[50];
  double R[50], T[50], res[50];
  double x[50], y[50], z[50];
  double px[50], py[50], pz[50], p[50];

  TFile *saveFile = new TFile(argv[1], "READ");
  TTree *saveTree = (TTree *)saveFile->Get("save");

  saveTree->SetBranchAddress("rawEvent", &event);
  saveTree->SetBranchAddress("nSeeds", &nSeeds);
  saveTree->SetBranchAddress("nSeedsX", &nSeedsX);
  saveTree->SetBranchAddress("nSeedsY", &nSeedsY);
  saveTree->SetBranchAddress("ax", ax);
  saveTree->SetBranchAddress("bx", bx);
  saveTree->SetBranchAddress("ay", ay);
  saveTree->SetBranchAddress("by", by);
  saveTree->SetBranchAddress("xIndex", xIndex);
  saveTree->SetBranchAddress("yIndex", yIndex);

  saveTree->SetBranchAddress("nCandidates", &nCandidates);
  saveTree->SetBranchAddress("time", &time);
  saveTree->SetBranchAddress("chisq", &chisq);
  saveTree->SetBranchAddress("charge", &charge);
  saveTree->SetBranchAddress("nRecHits", &nRecHits);
  saveTree->SetBranchAddress("nHodoHits", nHodoHits);
  saveTree->SetBranchAddress("iRecHits", iRecHits);
  saveTree->SetBranchAddress("dID", detectorID);
  saveTree->SetBranchAddress("res", res);
  saveTree->SetBranchAddress("R", R);
  saveTree->SetBranchAddress("T", T);
  saveTree->SetBranchAddress("x", x);
  saveTree->SetBranchAddress("y", y);
  saveTree->SetBranchAddress("z", z);
  saveTree->SetBranchAddress("px", px);
  saveTree->SetBranchAddress("py", py);
  saveTree->SetBranchAddress("pz", pz);
  saveTree->SetBranchAddress("p", p);

  TH1D *hist[24];
  for(int i = 0; i < 24; i++)
    {
      hist[i] = new TH1D(p_geomSvc->getDetectorName(i+1).c_str(), p_geomSvc->getDetectorName(i+1).c_str(), 100, -.2, .2);
      if((i & 1) == 0) 
	{
	  hist[i]->SetLineColor(kRed);
	}
      else
	{
	  hist[i]->SetLineColor(kBlack);
	}
    }

  int nEvtMax = saveTree->GetEntries();
  for(int i = 0; i < nEvtMax; i++)
    {
      saveTree->GetEntry(i);
      if(nRecHits < 17) continue;
      if(chisq > 100.) continue;
      if(p[0] < 25.) continue;

      bool hodoMask = true;
      for(int j = 0; j < 3; j++)
	{
	  if(nHodoHits[j] < 1)
	    {
	      hodoMask = false;
	      break;
	    } 
	}

      if(!hodoMask) continue;

      for(int j = 0; j < nRecHits; j++)
	{
	  hist[detectorID[j] - 1]->Fill(res[j]);
	}

      event->clear();
    }

  int order[24] = {3, 4, 1, 2, 5, 6, 10, 9, 11, 12, 7, 8, 16, 15, 18, 17, 14, 13, 21, 22, 19, 20, 23, 24};

  TCanvas *c1 = new TCanvas();
  c1->Divide(3, 4);
  for(int i = 0; i < 12; i++)
    {
      c1->cd(i+1);
      hist[order[2*i]-1]->Draw();
      hist[order[2*i+1]-1]->Draw("same");
    }

  c1->Print(argv[2]);

  double par[24];
  char buffer[300];
  fstream fin;
  fin.open("alignment.txt", ios::in);
  for(int i = 0; i < 24; i++)
    {
      fin.getline(buffer, 100);
      istringstream stringBuf(buffer);
      stringBuf >> par[i];
    }

  fin.close();
  
  fstream fout;
  fout.open(argv[3], ios::out);
  double z_plane[24], res_plane[24];
  double dz_plane[24], resol_plane[24];
  for(int i = 0; i < 24; i++)
    {
      cout << par[i] << "  " <<  hist[i]->GetMean() << "   " << hist[i]->GetRMS() << endl;
      fout << par[i] - hist[i]->GetMean() << "   " << hist[i]->GetRMS() << endl;
    
      z_plane[i] = p_geomSvc->getPlanePosition(i+1);
      dz_plane[i] = 0.;
      res_plane[i] = par[i] - hist[i]->GetMean();
      resol_plane[i] = hist[i]->GetRMS();

      //cout << i << "  " << z_plane[i] << "  " << par[i] << "  " << res_plane[i] << endl;
    }
  fout.close();

  TGraphErrors gr(24, z_plane, res_plane, dz_plane, resol_plane);
  //TCanvas *c2 = new TCanvas();
  //c2->cd(); gr.Draw("A*"); 
  gr.Fit("pol1");

  //c2->Print("test2.eps");

  delete c1;
  //delete c2;
  for(int i = 0; i < 24; i++)
    {
      delete hist[i];
    }

  return 1;
}
